Python Environment tests

# Basic Python
import os
import math
import numpy as np
# Opening Geospatial Data
import geopandas as gpd
import rasterio as rio
import rioxarray as rxr
# Plotting
import hvplot.xarray
import hvplot.pandas
import holoviews as hv 
# Other
import earthaccess
from modules.emit_tools import emit_xarray
from modules.ewt_calc import calc_ewt

# Hard Coded File urls
emit_url = 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20230401T203751_2309114_002/EMIT_L2A_RFL_001_20230401T203751_2309114_002.nc'
ecostress_url = 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/ECO_L2T_LSTE.002/ECOv002_L2T_LSTE_26860_001_10SGD_20230401T203733_0710_01/ECOv002_L2T_LSTE_26860_001_10SGD_20230401T203733_0710_01_LST.tif'

roi = gpd.read_file("../data/dangermond_boundary.geojson")

# Sign into EDL and get token
token = earthaccess.login(persist=True).token['access_token']

# Download EMIT Scene
# Get requests https Session using Earthdata Login Info
fs = earthaccess.get_requests_https_session()
# Retrieve granule asset ID from URL (to maintain existing naming convention)
granule_asset_id = emit_url.split('/')[-1]
# Define Local Filepath
fp = f'../data/{granule_asset_id}'
# Download the Granule Asset if it doesn't exist
if not os.path.isfile(fp):
    with fs.get(emit_url,stream=True) as src:
        with open(fp,'wb') as dst:
            for chunk in src.iter_content(chunk_size=64*1024*1024):
                dst.write(chunk)

Test EMIT tools reading and ortho

# Open EMIT and Orthorectify
emit_ds = emit_xarray(fp, ortho=True)

Geoviews Tiles, Image, and Polygon Plotting

# Geoviews + Polygon Visualizations
emit_ds['reflectance'].data[emit_ds['reflectance'].data == -9999] = np.nan
size_opts = dict(frame_height=405, frame_width=720, fontscale=2)
map_opts = dict(geo=True, crs='EPSG:4326', alpha=0.7, xlabel='Longitude',ylabel='Latitude')
emit_ds.sel(wavelengths=850, method='nearest').hvplot.image(cmap='viridis',tiles='ESRI',**size_opts, **map_opts)*roi.hvplot(color='#d95f02',alpha=0.5, crs='EPSG:4326')

# Crop and Save
emit_cropped = emit_ds.rio.clip(roi.geometry.values,roi.crs, all_touched=True)
emit_cropped.to_netcdf(f'../data/{emit_cropped.granule_id}_dangermond.nc')

Calculate EWT Using ray

The figure below should have a region with much higher ewt in the bottom left corner of the ROI. This likely indicates presence of iceplant (holds a lot of water).

%%time
# Calculate EWT
cropped_fp = '../data/EMIT_L2A_RFL_001_20230401T203751_2309114_002_dangermond.nc'
ds_cwc = calc_ewt(cropped_fp, '../data/', ewt_detection_limit=1.5, return_cwc=True)

2024-10-28 13:45:26,994 INFO worker.py:1816 -- Started a local Ray instance.
c:\Users\ebolch\repos\VITALS\python\modules\ewt_calc.py:117: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  data_vars = {"cwc": (list(emit_ds.dims.keys())[0:2], cwc.squeeze())}

# Visualize EWT
ds_cwc.hvplot.image(cmap='viridis',tiles='ESRI',**size_opts, **map_opts)

Open a COG using Bearer Token and Visualize Reprojected Data

# Open COG with Token
with rio.Env(GDAL_HTTP_AUTH="BEARER", GDAL_HTTP_BEARER=token):
    eco_ds = rxr.open_rasterio(ecostress_url, mask_and_scale=True, chunks='auto').squeeze('band', drop=True)
# Reproject to WGS84 and plot
eco_ds.rio.reproject('EPSG:4326').hvplot.image(x='x',y='y',**size_opts, cmap='inferno',tiles='ESRI', xlabel='Longitude',ylabel='Latitude',title='ECOSTRESS LST (K)', crs='EPSG:4326')

Test Interactive Plotting

Be sure to check a few points and ensure the spectral figure works.

# Mask out Bad Bands
emit_ds['reflectance'].data[:,:,emit_ds['good_wavelengths'].data==0] = np.nan

# Select RGB Bands
emit_rgb = emit_ds.sel(wavelengths=[650, 560, 470], method='nearest')

# Define and Apply Gamma Adjustment
def gamma_adjust(rgb_ds, bright=0.2, white_background=False):
    array = rgb_ds.reflectance.data
    gamma = math.log(bright)/math.log(np.nanmean(array)) # Create exponent for gamma scaling - can be adjusted by changing 0.2 
    scaled = np.power(np.nan_to_num(array,nan=1),np.nan_to_num(gamma,nan=1)).clip(0,1) # Apply scaling and clip to 0-1 range
    if white_background == True:
        scaled = np.nan_to_num(scaled, nan = 1) # Assign NA's to 1 so they appear white in plots
    rgb_ds.reflectance.data = scaled
    return rgb_ds

emit_rgb = gamma_adjust(emit_rgb,white_background=True)

# Interactive Points Plotting
# Modified from https://github.com/auspatious/hyperspectral-notebooks/blob/main/03_EMIT_Interactive_Points.ipynb
POINT_LIMIT = 10
color_cycle = hv.Cycle('Category20')

# Create RGB Map
map = emit_rgb.hvplot.rgb(fontscale=1.5, xlabel='Longitude',ylabel='Latitude',frame_width=480, frame_height=480)

# Set up a holoviews points array to enable plotting of the clicked points
xmid = emit_ds.longitude.values[int(len(emit_ds.longitude) / 2)]
ymid = emit_ds.latitude.values[int(len(emit_ds.latitude) / 2)]

first_point = ([xmid], [ymid], [0])
points = hv.Points(first_point, vdims='id')
points_stream = hv.streams.PointDraw(
    data=points.columns(),
    source=points,
    drag=True,
    num_objects=POINT_LIMIT,
    styles={'fill_color': color_cycle.values[1:POINT_LIMIT+1], 'line_color': 'gray'}
)

posxy = hv.streams.PointerXY(source=map, x=xmid, y=ymid)
clickxy = hv.streams.Tap(source=map, x=xmid, y=ymid)

# Function to build spectral plot of clicked location to show on hover stream plot
def click_spectra(data):
    coordinates = []
    if data is None or not any(len(d) for d in data.values()):
        coordinates.append(clicked_points[0][0], clicked_points[1][0])
    else:
        coordinates = [c for c in zip(data['x'], data['y'])]
    
    plots = []
    for i, coords in enumerate(coordinates):
        x, y = coords
        data = emit_ds.sel(longitude=x, latitude=y, method="nearest")
        plots.append(
            data.hvplot.line(
                y="reflectance",
                x="wavelengths",
                color=color_cycle,
                label=f"{i}"
            )
        )
        points_stream.data["id"][i] = i
    return hv.Overlay(plots)

def hover_spectra(x,y):
    return emit_ds.sel(longitude=x,latitude=y,method='nearest').hvplot.line(y='reflectance',x='wavelengths',
                                                                            color='black', frame_width=400)
# Define the Dynamic Maps
click_dmap = hv.DynamicMap(click_spectra, streams=[points_stream])
hover_dmap = hv.DynamicMap(hover_spectra, streams=[posxy])
# Plot the Map and Dynamic Map side by side
hv.Layout(hover_dmap*click_dmap + map * points).cols(2).opts(
    hv.opts.Points(active_tools=['point_draw'], size=10, tools=['hover'], color='white', line_color='gray'),
    hv.opts.Overlay(show_legend=False, show_title=False, fontscale=1.5, frame_height=480)
)

C:\Users\ebolch\AppData\Local\Temp\1\ipykernel_12556\166652677.py:8: RuntimeWarning: invalid value encountered in power
  scaled = np.power(np.nan_to_num(array,nan=1),np.nan_to_num(gamma,nan=1)).clip(0,1) # Apply scaling and clip to 0-1 range